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ABSTRACT 

■'  \ 

In  previous  work  by  the  author,  a  generalization  of  Godunov’s  method  for  systems  of  conservation  laws 
has  been  developed  and  analyzed  that  can  be  applied  with  arbitrary  time  steps  on  arbitrary  grids  in  one  space 
dimension.  Stability  for  arbitrary  time  steps  is  achieved  by  allowing  waves  to  propagate  through  more  than  one 
mesh  cell  in  a  time  step.  In  this  paper  the  method  is  extended  to  second  order  accuracy  and  to  a  finite  volume 
method  in  two  space  dimensions.  This  latter  method  is  based  on  solving  one  dimensional  normal  and  tangential 
Riemann  problems  at  cell  interfiices  and  again  propagating  waves  through  one  or  more  mesh  cells.  By  avoiding 
the  usual  time  stq>  restriction  of  explicit  methods,  it  is  possible  to  use  reasonable  time  steps  on  irregular  grids 
where  the  minimum  ceU  area  is  much  smaller  than  the  average  cell. 

Boundary  conditions  for  the  Euler  equations  are  discussed  and  special  attention  is  given  to  the  case  of  a 
Cartesian  grid  cut  by  an  irregular  boundary.  In  this  case  small  grid  cells  arise  only  near  the  boundary,  and  it  is 
desirable  to  use  a  time  step  appropriate  for  the  regular  interior  cells.  Numerical  results  in  two  dimensions  show 
that  this  can  be  achieved. 
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1.  Introduction. 

Hyperbolic  systems  of  conservation  laws  in  two  space  dimensions,  such  as  the  Euler  equations  of  gas 
dynamics,  have  the  form 

=0-  (11) 

Here  «  is  a  vector  of  m  variables  and  /  and  g  are  the  flux  functions.  In  order  to  solve  these  equations 
numerically,  the  domain  is  subdivided  into  computational  grid  cells  (or  intervals  in  one  dimension)  and  an 
ai^roximation  to  the  solution  over  each  cell  is  updated  in  each  time  step.  With  traditional  explicit  methods,  the 
appropriate  update  is  based  on  cell  values  in  only  a  few  neighboring  cells.  For  example,  with  a  three-point 
method  in  one  dimension,  only  the  nearest  neighbor  cells  on  either  side  can  effea  the  solution.  Most  finite 
volume  methods  in  two  space  dimensions  also  limit  the  domain  of  dependence  to  the  nearest  neighbors. 

This  imposes  an  obvious  limitation  on  the  time  step  due  to  the  CFL  condition,  which  states  that  a 
necessary  condition  for  convergence  is  that  the  domain  of  dependence  of  the  numerical  method  include  the 
domain  of  dependence  of  the  true  solution,  at  least  in  the  limit  as  the  grid  is  refined.  Fbr  a  three-point  method  in 
one  dimension,  this  means  that  the  time  step  k  must  be  restricted  by 

(1.2) 


where  is  length  of  the  smallest  cell  and  Smu  is  the  maximum  wave  speed  arising  in  the  problem  at 
hand.  The  quantity  v  defined  by  (1.2)  is  called  the  Courant  number. 

This  time  step  restriction  may  be  reasonable  on  a  unifonn  gri±  .\n  analysis  of  the  truncation  eirar  will 
often  show  that  (1.2)  is  consistent  with  the  choice  of  time  step  required  for  the  temporal  accuracy  to  agree  with 
the  spatial  accuracy.  On  the  other  hand,  if  the  grid  is  not  uniform  we  would  like  to  replace  A  gu  in  (1.2)  by  h  .y,, 
some  average  cell  length.  If  <k  A,,,  then  (1.2)  imposes  a  severe  time  step  restriction  that  may  greatly 
reduce  the  efficiency  of  the  method. 


This  situation  arises,  for  example,  in  certain  shock  tracking  procedures  (e.g.  [16]).  In  addition  to  a 
uniform  grid,  moving  points  are  introduced  that  may  fall  arbitrarily  close  to  the  fixed  points,  creating  small 
cells.  Similar  difficulties  are  encountered  in  die  moving  mesh  method  of  Harten  and  Hyman[12],  where  they 
must  take  special  care  to  insure  the  pr  i>er  separation  of  mesh  points. 


In  more  space  dimensions,  the  pr.  .^'^m  of  small  mesh  cells  is  difficult  to  avoid  in  several  conunon 
situations.  In  particular,  irregular  boundaries  cutting  through  a  regular  grid  (as  in  lugures  9.1  and  9.3)  result  in 
cells  with  arbitrarily  small  area.  We  would  like  to  choose  the  time  step  based  on  the  tegular  grid  cells  rather 
than  being  artificially  restricted  by  the  boundary  cells.  Small  cells  are  often  avoided  by  using  body  fitted 
coordinates,  so  that  the  grid  cells  have  more  uniform  areas.  However,  this  creates  additional  difficulties  in  first 
generating  and  then  computing  with  the  irregular  grids.  For  complicated  multicomponent  geometries, 
particularly  in  three  dimensions,  the  generation  of  such  grids  may  be  very  difficult,  requiring  complicated  grid 
patching.  This  has  lead  several  workers  to  reconsider  the  use  of  Cartesian  grids  (e.g.  [6,25,28]). 


<4 


Similar  difficulties  arise  if  a  shock  is  tracked  through  a  fixed  grid,  at  the  boundary  of  a  region  of  mesh  , 
refinement,  ch-  where  independent  grids  are  patched  together.  On  general  unstructured  grids,  such  as^  codi^ 
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In  this  paper  a  finite  volume  method  is  introduced  that  avoids  these  difficulties  by  dynamically  extending 
the  domain  of  dependence  in  such  a  way  that  the  CFL  condition  is  always  satisfied.  This  is  accomplished  by 
fuUy  exploiting  the  wave  propagation  that  is  inherent  in  hyperbolic  equations.  A  one  dimensional  Riemann 
problem  is  solved  normal  to  each  cell  interface  (meaning  the  initial  data  is  piecewise  constant  with  a  jump 
discontinuity  at  the  interface).  Solving  Riemann  problems  forms  the  basis  of  Godunov’s  classical  method  and 
numerous  finite  volume  methods,  but  in  most  such  methods  the  Riemann  solution  is  used  only  to  obtain  a  flux 
across  the  boundary.  This  flux  ic  then  used  to  update  only  the  neighboring  cell  values. 

In  the  method  proposed  here,  the  structure  of  the  Riemann  solution  is  used  to  a  much  greater  extent  The 
Riemann  solution  consists  of  waves  (shocks,  contact  discontinuities,  or  rarefaction  waves)  that  propagate  away 
from  the  interface,  each  carrying  an  increment  in  dte  flow  variables.  As  a  wave  propagates,  it  modifies  the  flow 
variables  in  the  cells  it  passes  through.  By  allowing  it  to  pass  through  more  than  one  cell,  the  CFL  conditicm 
remains  satisfied  with  larger  time  steps.  This  gives  a  generalization  of  Goduiwv’s  method  applicable  for  larger 
time  steps  on  arbitrary  grids.  In  one  space  dimension,  this  first  order  accurate  method  has  been  developed  and 
analyzed  in  a  series  of  papers  by  the  author[14,lS,161.  This  method  is  similar  in  spirit  to  the  "transport  collapse* 
method  of  Brenier[l,2,3].  In  the  present  paper,  the  method  is  extended  to  second  order  accuracy  and  arbitrary 
grids  in  two  space  dimensions. 

In  addition  to  propagating  waves  normally  from  cell  interfaces,  it  is  also  possible  to  incorporate 
tangential  wave  propagation  by  solving  auxilliary  Riemarm  problems  in  the  orthogonal  direction.  This  leads  to 
increased  stability  and  enhanced  modeling  of  two  dimensional  phenomena  and  propagation  of  shock  waves  that 
are  not  aligned  with  the  grid. 

An  additional  advantage  of  the  method  proposed  here  is  that  boundary  conditions  can  be  incoqx^rated 
into  the  algorithm  relatively  easily.  Farfield  nonreflecting  boundary  conditions  and  solid  wall  (or  moving 
piston)  boundary  conditions  for  the  Euler  equations  will  be  discussed  here. 

The  basic  idea  of  the  method  is  described  first  in  Section  2  for  a  scalar  equation  in  one  space  dimension. 
A  new  formulation  of  the  Lax-Wendroff  method  is  presented  and  then  converted  to  a  "high  resolution"  method 
(one  with  sharp  nonoscillatory  discontinuities  and  good  nonlinear  stability  properties).  The  formulation  of  this 
method  easily  allows  extension  to  arbitrary  grids  and  time  steps  and  iiKorporation  of  boundary  conditions. 

In  Section  3  the  method  is  extended  to  systems  of  conservation  laws  in  one  dimension.  Boundary 
conditions  for  the  Euler  equations  are  discussed  in  Section  4  and  then  some  one-dimesional  numerical  results 
are  presented  in  Section  S. 

In  the  remainder  of  the  piqrer  the  method  is  extended  to  a  finite  volume  method  in  two  dimensions.  The 
special  case  of  a  Cartesian  grid  with  irregular  boundaries  receives  particular  attention.  In  Section  9,  numerical 
results  are  presented  for  a  shock  tube  problem  oblique  to  the  grid  and  single  Mach  reflection  from  a  ramp. 

2.  One  dimensional  scalar  equations. 

We  first  consider  the  one-dimensional  scalar  conservation  law 

«*.+/(«).  =0.  (2.1) 

with  f{u)  assumed  to  be  a  convex  function.  The  computational  grid  is  given  by  xo<Xi  <  -  -  -  <Xm.  The 
mesh  spacing  can  be  arbitrary  and  we  set  A,-  =x,»  -x,-.  The  approximate  solution  on  the  mesh  cell  (x,  x,4.i) 
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time  is  denoted  by  t/*.  The  time  steps  can  also  be  arbitrary  and  varying.  Since  the  methods  are  one-step 
methods  and  will  be  described  over  a  single  time  step,  we  will  denote  the  current  time  step  simply  by 
k  =  -  f,  and  use  the  short  hand  abbreviations  £/;  ■  t/"  and  C/,  ■  t/"'*’*. 

For  each  value  of  j ,  set 


j(f  (Ujh)  -  Uj)  if  Uj  * 

Then  Sj  is  the  speed  at  which  the  discontinuity  between  Uj  and  I/y+i  propagates  according  to  the  Rankine- 
Hugoniot  relation.  For  convenience  at  this  point,  suppose  that  >  0  for  all  y .  Then  we  set 

\Uj  iatxKXj 


Uj(x  Jn)  =  '  Uj  +  (x- Xj^y^j  for  Xj^x< (2.3a) 
Uj^i  for  X  Sx^+i 

where  Xy+,/j=  -^(xy  +Xy+i)  and  Oy  is  the  slope  of  Uy  over  the  ;  th  interval,  as  yet  unspecified  (see  Figure  2.1). 

Based  on  these  functions,  we  define  an  approximate  solution  u(xj)  over  the  time  interval  ^  r  ^  At  time 
f,,  u(xj^)  is  simply  the  piecewise  linear  function  that  agrees  with  Uj(x,u)  on  the  ;  th  cell: 

u(x,/,)  =  «y(x,/,)  for  Xy  S X  < Xy^,  (2.4a) 

For  t  >  f,  we  let  the  wave  form  (2.3a)  propagate  with  speed  Sj,  defining 

Uy(x  J )  =  Hy(x  -  Sj(l-tJ,  u).  (2.3b) 

These  waves  are  then  combined  via  linear  superposition  to  obtain  the  approximate  solution 

u (x  y )  =  ^(x ,f,)  +  2  (.ujixj)  -  Uj{x ,tj).  (2.4b) 

i 

The  sum  in  (2.4b)  is  over  all  J  but  at  any  given  point  only  a  finite  number  of  terms  will  be  nonzero  since 
uy(x,r,)  is  constant  outside  of  (xy,zy4.,). 

The  s^rproximation  t/,-  at  time  is  obtained  by  projecting  the  approximate  solution  u(x  y,+i)  onto  the 
grid,  as  in  Godunov’s  method  (i.e.  by  averaging  over  the  mesh  cell  (x,  ,z,>i)): 


=  f/i  +  -rZr  (u/(xJ^^i)-Uj(xjJ)dx.  (2.5) 

rt.  j  * 

Note  that  U^  is  obtained  from  (/,  by  adding  to  it  contributions  from  several  neighboring  cells.  Figure  2.1  shows 
the  typical  contribution  to  indicated  by  the  shaded  area. 

For  a  linear  problem  u,  +  su,  =  0  with  coistant  speed  s,  the  approximation  u(xj)  is  in  fact  the  true 
solution  with  initial  data  u(x/«).  For  nonlinear  problems  the  waveform  (2.3a)  should  be  distorted  for  t  >  r,, 
and  so  (23b)  represents  a  local  linearization  of  the  flux  function.  The  superposition  (2.4b)  is  a  further 
linearization  of  the  nonlinear  problem.  In  spite  of  this,  it  is  possible  to  achieve  second  order  accuracy  with 
appropriate  choices  of  (Ty. 
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There  are  various  ways  to  specify  the  slope  Oj 

i.  Three  possibilities  are: 

II 

o 

(Godunov) 

(2.6a) 

-  Vj)  /  (Xy+V2  ~-*y+l/2) 

(Lax-WendrofO 

(2.6b) 

=  minmod  (af^,  a/f,) 

(minmod) 

(2.6c) 

where  the  minmod  function  is  defined  by 

ininmod(d^)=  ^(sgn(a)  +  sgn(6))  min(lal,ldf). 


The  names  assigned  to  the  first  two  choices  of  a  are  based  on  the  fact  that  the  method  reduces  to  one  of 
these  standard  methods  on  a  uniform  grid  with  Courant  number  v  <  1.  Then  each  is  influenced  by  at  most 
two  waves  and  the  method  can  be  rewritten  as  a  three-point  scheme.  Setting  v,  =  r./fe/A  and  computing  the 
integrals  in  (2.5)  explicitly,  we  obtain  (still  assuming  Sj>0  for  all  j) 

Vi  =  £/.  -  v.-,(f/i  -  £/._i)  -  v.-,)c,.,A  -  |vi(l  -  v.Oc.A.  (2.7) 

This  can  be  rewritten  in  conservation  form  as 

Ui^Ui-j^iFi-Fi.O  (2.8) 

where 

Fj=f(Uj)+^(l-Vj)sjajh.  (2.9) 

If  a  0  this  gives  the  first  order  upwind  scheme,  which  agrees  with  Godunov’s  method  for  the  problem 
considered  here.  With  slope  (2.6b),  this  is  easily  seen  to  be  the  standard  Lax-WendroCf  method. 

The  Lax-Wendroff  method  is  second  order  accurate  on  smooth  solutions,  but  it  generally  develops 
oscillations  and  possibly  instabilities  near  discontinuities  or  sharp  gradients.  The  minmod  choice  of  slopes 
(2.6c)  is  designed  to  avoid  this  difficulty.  Van  Leer[27]  used  a  similar  geometric  approach  with  limited  slopes 
to  define  his  second  order  accurate  MUSCL  schemes.  Later  Goodman  and  LeVeque[10,181  used  the  minmod 
slope  (2.6c)  to  derive  a  total  variation  diminishing  (TVD)  version  of  van  Leer’s  method.  The  method  (2.7)  with 
this  choice  of  a  is  closely  related  to  these  methods  and  to  a  variety  of  other  "high  resolution"  schemes  that  have 
recently  been  proposed  (e.g.  [1 1,19,20,29]).  In  particular,  for  a  linear  problem  with  f(u)  =  au  and  a  >  0,  (2.7) 
with  minmod  slopes  is  precisely  the  minmod  flux-limiter  method  introduced  by  Roe  (see  Sweby[26])  and  hence 
is  TVD. 

It  is  interesting  to  note  that  the  minmod  method  proposed  here  is  not  TVD  on  nonlinear  problems.  It  is 
possible  to  produce  examples  where  the  variation  increases  slightly,  although  detrimental  oscillations  or 
instabilities  have  never  been  observed.  Indeed,  the  lack  of  TVD  may  be  an  advantage  since  any  TVD  method 
must  reduce  to  first  order  accuracy  at  extreme  points[19].  In  numerical  experiments,  the  present  method 
appears  to  perform  better  at  such  points. 

The  main  advantage  ai  this  approach,  however,  is  that  we  can  easily  go  beyond  formula  (2.7)  and  apply 
the  method  as  described  by  (2S)  with  arbitrary  grids  and  time  steps. 
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Before  discussing  systems  of  conservation  laws,  we  point  out  a  slightly  different  formulation  of  the 
method  that  helps  simplify  its  implementation  and  extension  to  two  space  dimensions.  We  can  decouple  the 
wave  propagation  into  a  two  stage  procedure.  In  the  first  stage  we  propagate  the  wave  given  by  (2.3) 

with  Oy  =  0,  corresponding  to  the  first  order  algorithm.  In  the  second  stage  we  propagate  a  second  order 
correction  wave, 

«  ®(x  ,t )  =  Jiy(x  j)-up\xj) 

with  the  same  speed  Sj.  This  wave  has  integral  zero  and  initially  has  the  form 


uP\x.u)  =  \ 


0 


for  Xj^x  <  Xy+I 
otherwise 


When  viewed  in  this  way  it  is  clear  how  the  second  order  correction  relates  to  the  flux  corrected  transport 
(FCT)  algorithm  of  Boris  and  Book[3]  and  Zalesak[30]  and  to  the  ideas  of  Roe{22].  The  results  of  the  first 
order  accurate  scheme  are  corrected  by  redistributing  some  of  the  flux  between  grid  cells. 


In  the  above  discussion  we  assumed  that  Sj>0.  If  Sj  <  0  we  must  replace  (2.3)  by 


«>(x^)  = 


Vj 


Uj+(X-Xj^yj)aj 


for  X  <  Xj^i 
for  Xy+iSx  <Xj^2 
for  X  Sxy+2 


(2.10) 


in  order  to  obtain  a  method  that  reduces  to  the  Lax-Wendroff  or  minmod  method  in  the  special  cases  described 
above.  In  the  latter  case  we  also  redefine 


=  minmod  (a^’*’,  ojJ[). 

Note  that  the  linear  segment  in  Uy(x,r)  is  always  behind  the  discontinuity  as  it  propagates. 


3.  Systems  of  conservation  laws  in  one  dimension. 

We  now  develop  the  method  more  generally  for  hyperbolic  systems  of  conservation  laws  (2.1)  with 
u  €  R**.  We  assume  that  the  Jacobian  matrix  /'(«)  has  real  eigenvalues  X^(u)  for  all  u  and  that  each 
characteristic  field  is  either  genuinely  nonlinear  or  linearly  degenerate.  Then  the  Riemann  problem,  consisting 
of  (2.1)  with  piecewise  constant  initial  conditions 


X  <x,>, 
X  SXy*, 


(3.1) 


has  a  solution  for  -UjH  sufficiently  small.  This  similarity  solution  consists  of  m+1  constant  states 
separated  by  shocks,  contact  discontinuities,  or  rarefaction  waves.  (See  Lax[13]  for  a  description  of  the 
theory.) 


Although  the  method  can  be  adapted  to  use  this  exact  Riemann  solution,  we  will  base  the  method  on  an 
approximate  Riemann  solution  of  the  type  introduced  by  Roe[21].  The  Riemann  solution  of  (2.1)  with  data 
(3.1)  is  replaced  by  the  solution  of  the  linear  problem 


«,  +  AjU^  =  0 


(3.2) 


i 


where  Aj  ~Aj{Uj,Uj^i)  is  an  appropriate  average  value  of  /'(u).  This  gives  a  good  approximation, 
particularly  for  smooth  solutions  where  /£/ y+j  -  Ujll  =  O  (h ).  Denote  the  solution  of  this  problem  by  uj-'Kx  j ) 
for  t  S  Then  ,/)  consists  of  m  +  1  constant  states  separated  by  discontinuities  travelling  with  speeds 

JyiSjyjS  •••  (3.3) 

which  are  eigenvalues  of  Ay .  The  jump  in  across  the  p  th  such  wave  is  an  eigenvector  Rjp  of  A  j,  so  that 


and  £/y+i  -  f/y  is  decomposed  as 


-  ^Jp^Jp 


p“i 


(3.4) 


(3.5) 


To  obtain  high  resolution  we  will  also  introduce  a  correction  wave  up\x,t)  by  generalizing  the  scalar 
case  discussed  in  Section  2.  We  break  uf^  up  into  m  correction  waves  u-^.p  =  1,2,  ■■■  ,m  each  of  which 
travels  with  speed  and  has  an  initial  form  that  is  piecewise  linear  and  nonzero  over  a  single  mesh  interval. 
The  choice  of  interval  again  depends  on  the  sign  of  Sj^,: 


For  Sjp>Qr. 


r(jc-xy+i/2)<y> 

10 


for  Xy  S  X  <  Xy+i 

otherwise 


(3.6a) 


Fors^  SO: 


•(X 

0 


for  Xy*,  SX  <Xj^i 

otherwise 


tvhere  iiow  a  slope  vecu,’f  with  m  components,  tc  be  defined  below.  We  now  set 


i  uS^ixj)=  £  u^(x  -5^(t  -f,).tj 


and 


(3.6b) 


(3.7) 


My(X,r)=  b/*\x,/)  +  Uj^\xj). 


(3.8) 


Based  on  these  functions  Ujfxj),  we  can  define  an  approximate  solution  t((x,r)  just  as  in  the  scalar  case  via 
(2.4).  Averaging  u  (x  over  the  i  th  interval  as  in  (2.5)  gives  the  numerical  solution  (/,  at  time  r,.^, . 

If  (T^aO  then  Uj^\xj)»0  for  all  j.  In  this  case  the  method  reduces  to  the  large  time  step 
generalization  of  Godunov’s  method  (with  Roe’s  ^proximate  Riemann  solver)  tha  has  been  previously 
studied([14,lS,16]).  This  method  is  conservative  for  any  mesh  ratio  as  is  shown  in  [14]  where  the  method  is 
rewritten  in  conservation  form.  Note  that  the  method  remains  conservative  with  abitrary  slope  vectors  since 


f"  uy®(x,r)  dx  =  0  for  any  a^.  We  can  view  the  function  Uj{xj)  for  r  ^  r„  as  a  new  approximate  Riemann 

solver  which  gives  second  order  accuracy  f<a  an  appropriate  choice  of  a^.  Because  of  this  the  high  resolution 
method  with  nonzero  Ojp  can  also  be  written  in  conservation  form  following  Section  4  of  [14]. 

Corresponding  to  the  three  slope  choices  (2.6),  we  now  have 


(Godunov-Roe), 


(3.9a) 
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'  (^J*y2  -  Xj*i/7)  (Lax-WendrofO.  (3.9b) 

=  minmod  a^"')  (minmod),  (3.9c) 

where  in  (3.9c)  J  -  J  -  sgn(j^)  and  the  minmod  function  is  applied  to  vectors  componentwise.  Either  of  the 
last  two  choices  i^peais  to  give  a  second  order  accuracy  on  smooth  solutions.  In  order  to  achieve  stability  and 
nonoscillatory  behavior  near  discontinuities,  the  minmod  slope  (3.9c)  is  recommended  and  has  been  used  in  the 
numerical  experiments  reported  below. 

This  method  appears  difficult  to  implement  from  the  description  just  given.  In  fact,  it  is  quite  simple  if 
we  consider  each  wave  arising  &om  each  Riemann  problem  in  turn  and  compute  its  effect  on  the  solution.  Due 
to  the  linearization  of  wave  interactions  that  is  used  here,  each  wave  can  be  propagated  and  averaged  onto  the 
grid  independently  of  all  other  waves.  For  the  case  this  algorithm  has  already  been  described  in  some 

detail  in  [IS]  and  [16].  This  gives  propagation  of  the  waves 

The  correction  waves  simply  give  a  redistnbution  between  mesh  cells.  In  particular,  if  the  mesh  is 
uniform,  then  the  nonzero  portion  of  overlaps  at  most  two  cells,  independent  of  the  Courant  number,  and 
the  second  order  correction  is  easily  implemented. 

For  linear  systems  of  equations,  solution  of  the  Riemann  problems  gives  resolution  into  the  characteristic 
variables  and  u(xj)  is  again  the  exact  solution  with  the  piecewise  linear  data  u(z,/,),  for  all  r 
Consequently,  the  method  gives  second  order  accinacy  with  any  size  time  step  in  this  case.  For  nonlinear 
problems  with  smooth  solutions,  the  superposition  (2.4)  can  be  shown  to  be  consistent  with  second  order 
accuracy.  Towards  this  end,  a  semidiscrete  version  of  this  superposition  (discretized  only  in  time)  is 
investigated  in  [17]  and  shown  to  produce  an  approximate  solution  with  0(k^  error  over  a  time  step  of  length 
*. 


4.  Boundary  conditions  for  the  Euler  equations. 

One  advantage  of  the  approach  used  here  is  die  ease  with  which  boundary  conditions  can  be  derived  and 
implemented,  based  on  the  physics  of  the  problem.  This  will  be  illustrated  for  the  Euler  equations  at  a  free  field 
(non-reflecting)  boundary  and  at  a  solid  wall  {or  moving  piston)  boundary. 

The  Euler  equations  of  gas  dynamics  have  the  form 


+ 


(4.1) 


where  p,  v.p,  and  E  represent  the  density,  velocity,  pressure  and  energy,  respectively.  These  quantities  are 
related  by  £  =  p/(y-  1)  +  with  y=  1.4  for  air. 


Non-reflecting  boundary  conditions  are  easfly  achieved  by  simply  allowing  all  waves  which  cross  the 
boundary  to  pass  out  of  the  region,  with  no  incoming  waves.  In  a  subsonic  problem  these  boundary  conditions 
are  not  strictly  correct  For  example,  if  two  shocks  propagate  out  of  the  domain  at  different  times,  they  may 
interact  (in  the  exterior)  and  give  rise  to  a  weak  rarefaction  wave  which  should  reenter  the  domain  at  a  later 
time.  Correctly  modeling  such  behavior  requires  boundary  conditions  which  are  nonlocal  in  time.  Simply 
letting  the  waves  pass  out  appears  to  be  the  best  that  can  be  achieved  with  local  boundary  conditions  and  is 
easily  implemented  with  arbitrary  grids  and  time  steps. 
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Al  a  solid  wall  boundary  the  correct  boundary  condition  is  v  =  0.  This  gives  rise  to  a  boundary  Riemann 
problem  which  can  be  solved  for  the  incoming  wave.  Note  that  the  boundary  is  characteristic  in  this  case,  so 
that  there  is  only  one  incoming  wave. 

Another  way  to  impose  this  boundary  condition  for  the  PDE  is  to  extend  the  domain  beyond  the  boundary 
and  specify  initial  conditions  obtained  by  suitably  reflecting  the  values  from  the  interior.  For  example,  if  we 
wish  to  compute  over  time  k  and  represents  the  maximum  wave  speed,  then  it  suffices  to  extend  the 
domain  from  x  =  0,  say,  tox’  =  -ks^  and  set 

p(x)  =  p(-x) 

v(x)  =  -v(-x)  (4.2) 

p(x)  =  p(-x) 

for  X  *  ^  X  ^  0.  Solving  the  problem  over  this  extended  domain  (using,  for  example,  non-reflecting  boundary 
conditions  at  x  =  x*).  will  give  the  desired  solution  forx  S  0  with  the  boundary  condition  v(0,/)  =  0  satisfied. 

This  latter  approach  is  easy  to  implement  numerically  on  arbitrary  grids.  New  grid  points  are  introduced 
exterior  to  the  original  domain  by  reflecting  the  interior  points  which  fall  within  distance  ks^a  of  boundary. 
The  values  on  these  mesh  cells  are  obtained  by  reflecting  the  interior  values  as  in  (4.2).  Then  any  wave  leaving 
the  original  domain  during  the  time  step  is  matched  by  a  "reflected”  wave  which  enters  at  the  same  point  in 
time.  The  jump  in  v  in  the  entering  wave  is  the  negative  of  the  jump  in  v  in  the  exiting  wave  (due  to  the 
symmetry  in  solutions  to  the  Riemann  problem  under  the  reflection  (4.2))  and  hence  the  boundary  condition  is 
always  satisfied.  This  same  symmetry  exists  in  Roe’s  approximate  Riemann  solver  and  the  slopes  o^  will  also 
be  symmetric,  so  that  the  boundary  condition  remains  satisfied  for  the  method  discussed  here.  Again  the 
procedure  is  easily  implemented  for  arbitrary  time  steps  and  meshes.  Finally,  it  is  clear  that  the  method  has  the 
same  order  of  accuracy  near  the  boundary  as  in  the  interior,  since  after  the  reflection  the  boundary  essentially 
disappears  and  we  are  simply  using  the  interior  scheme  everywhere. 

If  the  solid  wall  boundary  is  moving  (as  for  example  in  a  piston  problem),  then  boundary  conditions  are 
also  easily  derived.  If  the  piston  is  located  at  x  =  z,  at  time  r,  and  is  moving  with  speed  So  (assumed  constant) 
for  ^  I  ^  then  the  mesh  cells  for  z,  ^x  ^  z,  -^s^k  are  reflected  to  the  left  of  z,  with  grid  values 

P(*ii  -  =  +Jt) 

-x)  =  2jo-v(z, +x)  (4.3) 

Pi^.-x)=p{x,+x) 

The  waves  generated  from  these  cell  boundaries  will  interact  with  the  interior  waves  in  such  a  way  that 
''(*•  +  -toff  ~  0  satisfied 

5.  Numerical  results  in  one  dimension. 

On  a  uniform  grid  with  v  <  1  the  present  minmod  method  performs  much  lilce  the  minmod  flux  limiter 
method  discussed  by  Sweby(261.  Figure  5.1  shows  some  results  for  a  shock  tube  problem  with  the  Euler 
equations  and  initial  data  with  v  =  0  everywhere  and  p  =  p  =  1  for  x  <  0.4,  p  =  p  =  3  for  x  >  0.4.  The  solution 
at  three  different  times  is  shown.  Note  that  the  shock  reflects  off  a  solid  wall  boundary  at  x  =  0  and  then 
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inieracis  with  the  contact  discontinuity.  Nonrefiecting  boundary  conditions  were  used  at  r:  =  1,  where  the 
rarefaction  wave  begins  to  leave  the  region. 

Figure  5.1a  shows  the  density  p  computed  using  the  first  order  accurate  method  with  Godunov  slopes 
=  0.  Figure  5.1b  shows  the  results  with  minmod  slopes.  In  each  case  h  =  0.02  and  k  =  h/2  corresponding 
to  V  =  .75.  The  solid  line  in  each  figure  is  the  solution  computed  with  201  mesh  points  using  minmod  slopes. 

Figure  5.2  demonstrates  the  ability  of  this  method  to  deal  with  nonuniform  meshes.  In  this  problem 
smooth  inidal  conditions  were  specified: 

p(x,0)  =  l.  +  0.1cos(2.5x),  v(x,0)  =  0.1sin(5x),  p(x,0)=  1.8  +  0.1cos(2.5x), 

again  with  a  solid  wall  boundary  at  x  =  0  and  nonreflecting  boundary  conditions  at  x  =  1.  The  solid  line  shows 
the  resulting  density  at  r  =  0.1  computed  with  51  equally  spaced  mesh  points  and  k  =  h  =  0.02,  giving  v  =  1. 
The  squares  show  the  solution  computed  with  51  mesh  points  Xj  =  (0  -1)^)^  densely  clustered  near  x  =  0. 
Again  k  =  0.02  but  now  the  Courant  number,  based  on  the  smallest  mesh  widtli  h ,,  is  roughly  2500.  The 
numerical  solution  does  not  appear  to  be  adversely  affected  by  this. 

6.  Finite  volume  methods  in  two  dimensions. 

In  two  space  dimensions,  the  wave  propagation  approach  can  be  used  to  rewrite  a  standard  finite  volume 
method  in  a  way  that  can  be  easily  extended  to  allow  larger  timesteps  than  are  normally  possible.  As  noted  in 
the  introduction,  this  might  be  valuable  with  nonuniform  grids  where  some  of  the  grid  cells  are  considerably 
smaller  than  average.  Standard  explicit  methods  have  a  stability  restriction  that  limits  the  time  step  in  terms  of 
the  smallest  cell  area.  Of  particular  interest  is  a  uniform  Cartesian  grid  in  which  some  of  the  grid  cells  are  cut 
off  by  an  irregular  boundary,  a  tracked  shock,  or  an  area  of  mesh  refinement  To  some  extend  the  discussion 
will  be  specialized  to  this  situation,  but  the  ideas  can  be  applied  on  more  general  unstructured  grids  as  well. 

In  one  space  dimension,  we  have  seen  that  nonuniform  grid  cells  are  easily  handled  by  allowing  waves  to 
propagate  entirely  through  a  small  cell  and  out  the  other  side.  The  same  approach  can  be  taken  in  two 
dimensions  with  waves  arising  at  cell  interfaces.  The  basic  algorithm  will  be  described  on  a  general 
unstructured  grid. 

Associated  with  each  grid  cell  Cy  is  an  approximation  Uj  to  the  solution  u{x^jJ  over  this  grid  cell. 
We  also  define  Ay  to  be  the  area  of  Cy  and  Ajy  to  be  the  length  of  the  face  separating  adjacent  cells  C,  and  C,. 
We  begin  by  describing  a  Goduiiov>type  scheme  and  will  later  introduce  high  resolution  corrections  based  on 
piecewise  linear  ^rproximations. 

In  a  finite  volume  method,  fluxes  across  the  cell  faces  are  defined  and  used  to  update  the  cell  values  on 
either  side  of  the  face.  Consider,  for  example,  the  face  separating  C,  and  Cy  in  Figure  6.1.  A  flux  can  be 
defined  by  solving  a  one-dimensional  Riemann  problem  normal  to  the  face.  Eiefine  new  variables  ^  (normal  to 
the  face)  and  t]  (tangential)  by 

^  =  ax  +  Py,  Ti  =  -Px  +  ay  (6.1) 

with  o  =  COS0, 3  =  sin0  where  0  is  the  angle  of  the  cell  interface.  The  conservation  laws 


«!+/(«),+«(“),  =0 


(6.2) 


-lo¬ 


an;  transformed  to 


u,  +F(u\  +  G(u\  =  0 

where 


(6.3) 


f(u)  =  a/(u)+P«(u)  (6.4) 

G(u)  =  -P/(tt)  +  a^(H). 

In  the  new  coordinates  u  is  constant  in  t]  and  so  (6.3)  reduces  to  a  one-dimensional  Riemann  problem 

u,+F(u)^  =  0  (6.5) 

with  left  and  right  states  (/,  and  Uj.  For  rotationally  invariant  equations  such  as  the  Euler  equations  this  is 
particularly  simple  since,  after  a  change  of  dependent  variables  to  rotate  the  velocity  field,  the  form  of  F  (u ) 
agrees  with  /  (u)  and  hence  a  single  Riemann  solver  suffices  for  all  angles  8. 

As  in  the  one-dimensional  algorithm,  we  use  Roe’s  approximate  Riemann  solver  and  obtain  a  set  of 
waves  traveling  with  speeds  Si^2‘' ' '  -^m  in  the  ^-direction.  As  before  we  denote  the  jumps  in  u  across  these 
waves  by  the  vectors R  i,  ••  •,/?,«,  so  that 


r 


(6.6) 


Figure  6.1  shows  these  waves  propagating  from  the  interface  in  a  situation  where  the  waves  lie  entirely  within 
neighboring  cells. 


We  think  of  u  being  incremented  by  -sgnfr^jR^  at  each  point  swept  out  by  the  ^th  wave.  Consequently, 
the  average  value  U,  of  u  over  the  cell  C,  in  Figure  6.1,  for  example,  will  be  incremented  by 


If jlkfiy 


(6.7) 


The  factor  multiplying  R  i  is  the  ratio  of  the  area  swept  out  by  the  1-wave  to  the  total  cell  area.  Similarly,  Uj  is 
incremented  by 


(-Rj)  + 


(6.8) 


These  increments  to  Uj  and  Uj  can  be  interpreted  more  traditionally  in  terms  of  a  flux-difference 
splitting.  The  total  increment  to  u  (which  has  been  split  between  Q  and  C,  and  averaged  over  these  cells)  is 
the  sum  of  A,-  times  (6.7)  plus  A ^  times  (6.8),  giving 

-Wl,y[fjRl  +  f2R2  +  ^3/?3].  (6.9) 

A  basic  property  that  Roe’s  approximate  Riemann  solver  shares  with  the  exact  Riemann  solution  is  that 


2f,R^=f(f/;)-F(£/.) 


so  that  (6.9)  becomes 


(6.10) 


-khijiF{Uj)-FiUi)). 


(6.11) 
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This  is  ihe  flux  difference  between  C;  and  Cy  which  has  been  split  via  (6.7)  and  (6.8). 

The  method  as  described  so  far  is  thus  equivalent  to  traditional  flux-difference  splitting  methods.  By 
integrating  the  conservation  laws  around  each  cell,  it  is  easily  shown  that  such  methods  are  conservative  for  any 
splitting  of  the  flux  difference.  A  variety  of  such  methods  have  been  proposed,  most  of  which  share  the 
property  that  flux  differences  are  split  only  between  and  Cj  with  no  contribution  to  other  cells.  This  leads  to 
a  restriction  on  the  time  step,  since  clearly  the  factor  multiplying  R ,  in  (6.7),  for  example,  should  be  no  greater 
than  1  for  stability.  If  A,-  is  very  small,  this  gives  a  severe  restriction  on  k . 

The  advantage  of  the  wave  propagation  interpretation  of  Figure  6.1  is  twofold.  Fust,  it  makes  it  clear  that 
in  some  circumstances  other  cells  should  also  be  affected.  More  importantly,  it  makes  it  easy  to  determine 
exactly  what  that  e^ect  should  be.  Figure  6.2  shows  a  typical  wave  propagating  over  a  different  grid.  Each  ceil 
value  should  clearly  be  updated  by  sgn(s,)l?,  times  the  ratio  of  the  shaded  area  within  the  ceil  to  the  total  cell 
area.  Note  that  this  ratio  is  always  less  than  or  equal  to  1.  Splitting  each  wave  in  this  manner  and  adding  up  all 
increments  to  all  cells  gives  the  same  total  increment  (6.9)  and  hence  the  method  remains  conservative,  while 
being  much  more  stable  in  a  situation  such  as  Figure  6.2  than  a  method  that  increments  only  Cj  and  Cy. 

Instabilities  can  still  arise,  however,  as  illustrated  by  the  following  simple  example.  Consider  the  scalar 
linear  advection  equation 

«,-m,  +  H,=0  (6.12) 


on  a  uniform  rectangular  mesh  with  initial  data  consisting  of  a  sawtooth  discontinuity  at  45”  to  the  mesh; 


ifi +;■ 


for  some  M .  U  k  ^hl2  then  the  method  described  above  is  stable  and  in  fact  fa  k  -  hl2  gives  the  exact 
solution,  a  wave  propagating  through  one  additional  diagonal  of  mesh  cells  in  each  time  step: 


_  Jl  iii+j^M  +  l 
ifi+jiM  +2 

For  larger  k,  however,  the  method  is  unstable.  For  example,  ilk  =  h  then  the  wave  should  propagate  through 
two  additional  diagonals  in  each  step.  Instead  the  waves  propagating  in  the  x-  and  y -directions  overlap  and 
give  the  value  Uij  =  2for  i  +J  =  M  +  l  with  no  propagation  into  the  second  diagonal.  In  later  time  steps  this 
instability  grows  exponentially. 


In  order  to  avoid  this  difficulty,  it  appears  to  be  necessary  to  include  tangential  wave  propagation  as  well 
as  normal  propagation.  For  equation  (6.12)  it  is  easy  to  see  how  to  include  tangential  propagation.  All  waves 
should  propagate  with  tangential  velocity  1  as  weU  as  normal  velocity  1.  Figure  6.3  shows  this  wave 
propagation  in  the  case  k  =h/2,  giving 


T 

0.75 

0.25 

0 


if  i+j<M  +  l 
if  i  +  y  =  M  +  1 
ifi  +  y  =  Af  +2 
if  i+j>M +  2 
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Thc  method  now  gives  exact  propagation  with  k  =  h  and  indeed,  except  for  smearing  introduced  by  the 
averaging  process,  this  modified  method  gives  exact  propagation  for  any  time  step  k . 

For  nonlinear  systems  of  equations  the  tangential  propagation  is  introduced  as  follows.  After  solving  the 
normal  Riemann  problem  we  have  waves  traveling  with  speeds  normal  to  the  cell  face.  Suppose  now  that 
we  split  up  the  wave  into  subwaves  R^,  •  •  • ,  R^’^^  with 

(6.13) 

«-« 

and  assign  tangential  velocities  •  •  • ,  The  wave  in  Figure  6.2,  for  example,  might  be  split  into  two 
waves  R^^^  and  R^  with  tangential  speeds  and  s,®  as  shown  in  Figure  6.4.  Note  that  the  area  swept  out 
by  each  wave  is  I  s,  I  khij,  the  same  as  the  area  of  the  original  wave  in  Figure  6.2,  indqiendent  of  the  tangential 
speeds  If  each  of  these  waves  is  propagated  as  before  to  increment  neighboring  mesh  cells,  the  total 
increment  will  be 

-kh,j[s^R<^^  +  s,/?  =  -khijS^Rp  (6.14) 

in  view  of  (6.13).  Hence  the  method  remains  conservative  for  an  arbitrary  splitting  R^^^  and  arbitrary  speeds 
provided  only  that  (6.13)  is  satisfied. 

A  natural  way  to  choose  the  splitting  is  by  solving  a  tangential  Riemann  problem  with  initial  data  having 
a  discontinuity  of  magnitude  Rp.  For  the  equation  we  take  the  tangential  portion  of  (6.3), 

«.+G(u)^  =  0.  (6.15) 

Figure  6.1  suggests  a  possible  set  of  initial  data.  The  wave  Rp  propagates  with  the  state 

Ui  +  L  (6-16a) 

i<p 

to  one  side  and 

Uj~^R^  (6.16b) 

n>p 

to  the  other.  Taking  these  as  the  left  and  right  states  for  the  equation  (6.15)  gives  a  discontinuity  of  magnitude 
Rp  which  will  then  be  resolved  into  waves  R^^^,  ,  R^"^  propagating  in  the  q  direction  with  speeds 

,0)  ...  r  (") 

*  f  • 

For  the  linear  equation  (6. 12),  m  =  1.  s  j  =  1  and  s  =  1  for  all  t , ; .  This  gives  the  correct  propagation 
as  described  above.  More  generally,  consider  a  linear  system 

u,  +Au„  +BUy  =  0. 

First  suppose  that  A  and  B  have  identical  eigenvectors  and  hence  are  simultaneously  diagonalizable.  Then 
each  wave  Rp  obtained  by  splitting  Uj  -  Ui  into  eigenvectors  of  A  is  already  an  eigenvector  of  S  as  well,  so 
that  R^'  =  Rp  and  R^^^  =  0  for  q  The  normal  and  tangential  speeds  Sp  and  are  the  corresponding 
eigenvalues  of  A  and  B  respectively.  In  this  case  the  method  gives  exact  propagation  (followed  by  smearing 
due  to  the  averaging  process)  for  any  time  step  k  and  we  are  essentially  using  the  method  of  characteristics. 
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If  A  and  B  arc  noi  simultaneously  diagonalizabte,  or  for  nonlinear  systems,  each  wave  will  in  general 
be  split  into  m  waves  Except  in  special  circumstances  the  method  will  not  give  exact  propagation  and  the 
time  step  will  need  to  be  restricted  in  order  to  obtain  reasonable  results.  Also  note  that  when  waves  from 
din^erent  interfaces  overlap,  they  are  combined  by  linear  superposition  as  in  the  one  dimensional  algorithm. 
This  incurs  additional  errors  on  nonlinear  problems.  A  serious  error  analysis  of  the  method  has  not  yet  been 
attempted,  but  the  numerical  results  of  Section  9  give  some  indication  of  its  capabilities. 

In  the  special  case  of  a  uniform  Cartesian  grid,  with  h^j  m  h  and  A,  s  h  the  wave  propagation  becomes 
relatively  simple,  particularly  if  the  Courant  number  is  less  than  1.  Then  each  wave  affects  at 

most  two  cells,  as  in  Figure  6.3,  and  the  areas  of  intersection  are  easy  to  compute.  Although  we  don’t  present 
the  details,  the  method  can  be  written  as  a  nine-point  finite  difference  method  in  conservation  form.  If 
tangential  wave  propagation  is  suppressed,  the  method  reduces  to  a  five-point  scheme. 


7.  High  resolution  in  two  dimensions. 

Greater  accuracy  for  smooth  flows  and  better  resolution  of  discontinuities  can  be  obtained  by  introducing 
piecewise  linear  approximations  as  in  the  one  dimensional  method.  Ideally  one  might  like  to  use  bilinear 
approximations  matching  slopes  in  the  x-  and  > -directions  simultaneously,  but  it  turns  out  that  a  simpler 
approach  is  sufficient  to  reproduce  the  Lax-Wendroff  method  for  scalar  problems  on  uniform  grids. 
Introducing  minmod  slope  limiting  gives  a  stable  second  order  correction  that  is  easy  to  generalize  and 
implement  for  nonlinear  systems  on  arbitrary  grids. 

First  consider  a  uniform  Cartesian  grid  with  Courant  number  less  than  one  and  suppress  tangential 
propagation.  Then  we  obtain  a  five-point  flux-splitting  method  as  discussed  above.  We  can  think  of 
implementing  this  by  applying  the  first  order  (Hie  dimensional  algorithm  of  Section  3  simultaneously  in  the  x- 
and  y -directions.  It  is  natural  then  to  also  apply  the  second  order  correction  waves  inutxluced  there.  In 
handling  the  interface  between  Uij  and  for  example,  we  obtain  waves  i,  ■  ■  ■  ,R^  propagating  in  the 
x-directi(Hi.  Suppose  that/?^  propagates  with  speed  >  0.  Then  is  incremented  by 


according  to  the  first  order  algorithm.  For  the  second  order  conectiem  we  propagate  the  wave  u  ^  fiom  Section 
3,  viewed  now  as  a  two  dimensional  wave  that  is  constant  in  y  over  the  grid  cells  in  question  and  piecewise 

linear  in  x.  Thus  Un.ij  is  incremented  by  -  ks^yCp  while  C/jy  is  incremented  by 

-kSp)af,  where  o,  is  the  slope  vector  in  the  x -direction,  again  defined  by  one  of  (3.9).  Wave 

2  H 

propagation  in  the  y  -direction  is  handled  in  an  analogous  manner. 

If  we  choose  the  Lax-Wendrofif  slope  (3.9b),  <T^  =Rplh,  then  for  scalar  equations  the  methcxl  just 
desexibed  reduces  to  the  following  two  dimensional  version  of  the  Lax-WendroH'  method: 


ai) 
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(£iUij*xyg(Vij)f  ^  (£(UijygiUij.x)f' 

■"  ■"  U,rUiJ-X 

This  method  is  second  order  accurate.  As  before,  in  order  to  eliminate  oscillations  it  is  better  to  use  the 
minmod  slope  (3.9c)  in  practice,  giving  a  method  that  remains  second  order  accurate  except  at  extrema. 

When  tangential  propagation  is  included,  it  is  possible  to  split  up  the  second  order  correction  wave  into 
subwaves  and  propagate  these  in  the  tangential  direction  as  well.  Here  we  adc^t  a  simpler  approach  and  always 
propagate  the  correction  wave  normal  to  the  interface.  On  uniform  grids  the  second  order  correction  is  then 
very  easy  to  implement  and  adds  little  cost  to  the  algorithm.  The  numerical  results  of  Section  9  indicate  that 
this  is  sufficient  to  give  increased  resolution  of  two-dimensional  phenomena. 

On  unstructured  grids,  a  high  resolution  correction  can  be  defined  in  various  ways.  Here  is  one  untested 
possibility,  motivated  by  the  above  treatment  of  uniform  grids.  Let  4  and  t)  be  the  local  coordinate  directions 
defined  by  the  interface  between  Q  and  Cy.  For  the  pth  wave,  define  a  slope  vector  in  the  ^-direction 
(normal  to  the  interface).  The  simplest  choice  would  be  the  "Lax-WendrofT  slope 

<5,=R,lh'  0-2) 

where  A*  is  some  measure  of  the  normal  distance  between  ceils  C,  and  Cj,  such  as  the  difference  in 
coordinates  of  the  centers  of  mass  of  these  two  cells.  Then  a  piecewise  linear  wave  with  this  slope  in  the 
direction,  zero  slope  in  the  r^-direction.  and  total  integral  zero  can  be  propagated  in  the  ^-direction  with  velocity 
s,.  The  shape  of  the  region  where  the  slope  is  is  still  undetermined.  Following  the  one  dimensional 
method,  we  could  choose  the  region  to  be  the  cell  Cj  if  s,  >  0  and  Cy  if  <  0.  To  be  specific,  assume  that 
r,  >  0  and  let  (^,r|f)  be  the  center  of  mass  of  the  ceil  C,-.  in  the  local  coordinates  determined  by  the  interface. 
Then  the  high  resolution  correction,  which  must  be  averaged  onto  the  cells,  is  the  piecewise  linear  function 

f(^  -  Spk  -  if  (^  -  5, A  ,Ti)  G  Ci 

“  (^.'H.^ii+i)  -  |q  otherwise. 

In  problems  with  shocks  one  would  again  like  to  introduce  a  slope  limiter  instead  of  using  (7.2)  directly. 
The  best  way  to  do  this  on  unstructured  grids  has  not  yet  been  determined. 

8.  Boundary  conditions  in  two  dimensions. 

Nonreflecting  outflow  boundary  conditions  are  implemented  in  two  dimensions  in  exactly  the  same  way 
as  in  one  dimension.  Waves  propagating  out  of  the  computational  domain  at  such  a  boundary  simply  disappear 
and  no  new  waves  are  propagated  inward  from  boundary  interfaces. 

At  a  solid  wall  boundary,  the  pnsper  boundary  conditions  for  the  Euler  equations  are  obtained  by 
specifying  zero  normal  velocity.  If  the  boundary  is  flat,  this  boundary  condition  can  be  implemented  in  the 
same  way  as  was  done  in  Section  4  for  the  one  dimensional  Euler  equations.  We  will  consider  this  simple  case 
first  and  see  how  to  reintennet  it  in  a  manner  more  suited  to  arbitrary  boundaries.  For  flat  boundaries,  the  grid 
cells  near  the  boundary  can  be  reflected  across  the  boundary  to  form  fictitious  cells  as  in  Figure  8.1.  We  assign 
values  of  C/  in  these  exterior  cells  based  on  the  values  in  the  mirror-image  interior  cells  so  that  the  density, 
pressure,  and  tangential  velocity  in  the  exterior  ceil  agree  with  those  values  in  the  interior  while  the  normal 
velocity  is  negated.  Then  we  can  simply  compute  over  the  extended  domain  using  the  algorithm  described 
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above,  ignoring  the  boundary.  Each  wave  from  the  interior  which  crosses  the  boundary  will  be  matched  by  a 
wave  from  the  exterior  in  such  a  way  that  the  normal  velocity  in  mirror-image  cells  adjacent  to  the  boundary 
will  remain  negatives  of  one  another,  simulating  zero  normal  velocity  on  the  wall  itself. 

Alternatively,  rather  than  creating  fictitious  exterior  cells  and  propagating  waves  from  their  interfaces  to 
match  interior  waves,  we  can  simply  reflect  the  interior  waves  when  they  hit  the  boundary  as  suggested  by 
Figure  8.1.  If  the  wave  carries  increments  Rp  then  the  reflected  wave  carries  the  same  increments  of  density, 
pressure,  and  tangential  velocity.  The  increment  of  normal  velocity  is  negated. 

The  latter  approach  has  the  advantage  that  it  extends  more  easily  to  general  boundaries.  If  the  boundary 
is  not  flat,  creating  globally  well-defined  fictitious  cells  by  reflecting  the  interior  cells  is  not  possible.  It  is, 
however,  still  possible  to  reflect  waves  at  the  boundary  provided  only  that  we  assume  the  boundary  is  locally 
flat  where  it  intersects  the  wave,  as  can  be  done  to  good  approximation  if  the  boundary  is  smooth. 

In  addition  to  handling  the  reflection  of  waves  arising  firom  interior  interfaces,  it  is  also  necessary  to 
generate  a  wave  from  the  cell  face  defining  the  boundary  itself.  We  solve  a  one  dimensional  Riemann  problem 
normal  to  the  boundary.  The  left  state  comes  from  the  interior  celL  The  right  state  is  the  same  except  that  the 
normal  velocity  is  negated,  as  motivated  by  Figure  8.1  and  the  previous  discussion.  Because  the  two  states  are 
related  in  this  manner,  the  solution  of  the  Riemann  problem  will  give  only  one  wave  propagating  into  the 
computational  domain,  with  zero  normal  veloci^  in  its  wake.  Of  the  other  three  waves,  two  will  have  zero 
velocity  and  the  third  will  propagate  out  of  the  computational  domain  and  is  ignmed. 

9.  Numerical  results. 

To  date  the  method  has  been  implemented  only  on  a  Cartesian  grid  with  "small"  time  steps,  meaning  the 
Courant  number  is  less  than  one  relative  to  the  area  A  ^  of  the  regular  grid  cells.  However,  boundaries  are 
allowed  to  cut  through  this  grid  and  the  area  of  the  resulting  boundary  cells  may  be  arbitrarily  small.  The 
boundary  is  approximated  by  piecewise  linear  segments  connecting  the  points  where  the  boundary  crosses  the 
grid  lines.  Consequently  aU  cells  are  polygons  with  at  most  five  sides,  since  we  assume  each  cell  contains  at 
most  one  boundary  segment 

Because  of  the  Courant  number  restriction,  wave  propagatioa  in  the  interior  cells  is  easily  implemented 
as  discussed  in  Section  6.  It  is  only  near  the  boundaries  that  waves  may  intersect  several  cells  in  an 
unpredictable  manner.  Each  wave  is  also  a  polygon  which  can  be  represented  by  a  list  of  its  vertices.  This 
polygon  can  be  intersected  with  the  polygonal  grid  cells  and  the  area  of  intersection  calculated  using  standard 
algorithms  of  computational  geometry.  A  wave  that  strikes  the  boundary  is  reflected  by  intersecting  it  with  a 
half  space  which  locally  approximates  the  exterior.  This  gives  a  polygon  representing  the  portion  of  the  wave 
exterior  to  the  domain  which  is  then  reflected  through  the  boundary  of  the  half  space  to  obtain  an  interior 
polygon. 

This  basic  algorithm  has  been  fully  implemented.  The  high  resolution  correction  terms  have  been 
properly  included  only  at  the  regular  grid  cells  in  the  interior.  At  the  boundary  an  approximation  is  used  that 
gives  reasonable  results  but  can  probably  be  improved  upon. 

The  test  cases  presented  here  involve  only  flat  boundaries.  The  handling  of  boundary  conditions  is  then 
equivalent  to  reflecting  the  interior  cells  across  the  boundary.  This  minimizes  the  possible  negative  eHiects  of 
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the  boundary  conditions  and  gives  a  reasonable  test  of  the  basic  algorithm  in  situations  where  there  are  irregular 
small  grid  cells. 

Allowing  arbitrarily  small  boundary  cells  poses  some  numerical  difficulties  not  associated  with  the 
theoretical  stability  properties  of  the  method.  If  the  size  of  a  grid  cell  is  on  the  order  of  the  machine  rounding 
error,  then  the  area  of  the  intersection  of  that  cell  with  a  wave  cannot  be  accurately  computed.  This  can  lead  to 
gross  errors  that  will  eventually  cause  noticeable  errors  in  other  grid  cells  and  possibly  nonphysical  flow 
quantities.  In  order  to  avoid  such  problems,  we  have  chosen  to  ignore  any  boundary  grid  cells  with  area  less 
than  10"^h^  This  amounts  to  introducing  imperceptible  holes  in  the  otherwise  solid  boundary.  The 
experiments  reported  here  were  performed  on  a  Vaxstation  II  in  single  precision. 

The  first  test  case  is  a  one  dimensional  shock  tube  problem  in  a  channel  at  40°  to  the  grid.  We  take  the 
same  initial  data  as  in  Section  S.  Figure  9.1  shows  the  boundary  cells  of  the  computational  grid  and  contour 
lines  of  constant  density  for  the  numerical  solution  at  t  =  0.3,  using  the  high  resolution  method  with  h  =  1/40 
and  timestep  k  =  0.007S.  This  gives  a  Couiant  number  of  roughly  0.S  relative  to  the  regular  cells,  though  it  is 
nearly  100  times  larger  for  the  smallest  irregular  cells  that  are  not  ignored.  Figure  9.1a  shows  the  results 
obtained  with  the  basic  first  order  method.  Figure  9.1b  shows  the  results  when  the  high  resolution  correction 
waves  discussed  in  Section  7  are  included.  Both  figures  show  contour  lines  of  density  as  well  as  the  irregular 
grid  cells  adjacent  to  the  boundary.  Hgure  9.2  shows  corresponding  plots  of  the  density  along  the  lower  wall  of 
the  channel.  Values  in  each  boundary  cell  are  plotted  vs.  distance  along  the  walL  No  attempt  has  been  made  to 
extrapolate  to  the  actual  boundary  location,  or  to  account  for  the  nonuniformity  of  these  boundary  cells. 
Consequendy  some  nonsmoothness  is  visible  in  diis  boundary  data  (and  even  more  so  in  Figure  9.4  below). 
Plots  of  grid  values  down  the  center  of  the  channel  would  show  much  smoother  results. 

Smearing  is  clearly  reduced  by  introducing  the  high  resolution  correction  waves.  In  both  cases  the 
contour  plots  show  that  the  method  is  capable  of  maintaining  one-dimensionality  of  the  solution,  with 
essentially  the  same  accuracy  in  the  irregular  cells  as  in  the  interior.  The  small  boundary  cells  cause  no  stability 
problems. 

The  second  test  example  is  an  unsteady  shock  reflection  from  a  27°  ramp  at  Mach  2.03,  with  single  Mach 
reflectioiL  Experimental  as  well  as  numerical  results  for  this  problem  can  be  found  in  the  paper  of  Glaz. 
CoUela,  Glass,  and  Deschambault  (Case  2  in  [9]).  The  irutial  flow  conditions  are  p  =  0.387,  v  =0,p  =  33.3 
ahead  of  the  shock  and  p  =  l.OS,  v  =  -14.06,p  =  1S4.S  behind  the  shock,  where  v  if  the  velocity  normal  to  the 
shock.  The  shock  was  initially  distance  0.05  firoro  the  ramp  comer  and  results  are  presented  at  time  0.27S  with 
h  =  1/80  and  k  =  2J5xl0~*  (1 10  time  steps).  The  results  presented  here  were  computed  with  the  high  resolution 
method.  Computations  with  the  first  order  version  have  also  been  done  and  give  the  same  qualitative  results 
with  much  more  smearing. 

For  comparison  purposes,  two  different  orientations  of  the  ramp  are  considered.  In  Case  I  (Hgure  9.3a), 
the  grid  is  aligned  with  the  ramp.  This  is  the  orientation  often  used  in  other  codes  to  avoid  irregular  grid  cells  in 
the  region  of  primary  interest  (e.g.  [9],  [29]).  In  Case  II  (Figure  9.3b)  the  grid  is  aligned  with  the  incident  shock 
aixl  the  ramp  cuts  through  the  grid  obliquely.  In  each  case  30  equally  spaced  density  contours  are  plotted.  The 
results  agree  closely  and  have  essentially  the  same  resolution. 

Bgure  9.4  shows  a  plot  of  relative  density  (p/po.  where  po  =  0.387)  along  the  lower  walL  The  solid  line 
shows  the  results  from  Case  I  while  the  +  symbols  show  results  from  Case  H  The  squares  are  experimental 
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resulls  taken  from  [9].  nie  results  are  in  good  agreement  and  also  compare  reasonably  well  with  numerical 
results  of  [9].  The  discontinuities  obtained  here  arc  not  as  sharp,  particularly  the  slip  line  connecting  the  triple 
point  to  the  wall.  In  large  part  this  is  due  to  the  coarseness  of  the  grid  (essentially  80x40  grid  points  compared 
with  350x130  in  [9]).  Figure  9.4  shows  that  even  in  Case  H,  where  the  ramp  is  oblique  to  the  grid,  there  arc  at 
most  two  points  in  the  shocks. 

10.  Conclusions. 

An  approach  for  deriving  high  resolution  numerical  methods  has  been  presented  that  allows 
generalization  to  stable  methods  on  arbitrary  grids.  Numerical  results  in  two  space  dimensions  show  that  the 
accuracy  in  irregular  grid  cells  is  comparable  to  what  is  obtained  with  regular  cells.  In  particular,  for  Car"‘sian 
grids  cut  by  inegular  boundaries,  it  is  possible  to  obtain  stable  high  resolution  results  with  a  time  step  chosen 
based  on  the  interior  regular  cells. 

Another  application  of  these  ideas  is  to  mesh  refinement  using  overlapping  rotated  Cartesian  grids,  as  in 
the  work  of  Berger(l,2].  Again  irregular  cells  are  created  at  the  grid  interfaces  and  the  present  approach  gives  a 
way  of  deriving  conservative  interface  conditions.  This  is  currently  being  studied  by  the  author  and  M.  Berger. 

A  similar  approach  toward  achieving  stability  has  been  taken  by  Chem  and  Colella[S].  They  give  an 
algorithm  for  shock  tracking  across  a  Cartesian  grid  in  which  fluxes  are  restributed  from  small  grid  cells  to 
neighboring  grid  ceils  to  insure  stability.  Their  redistribution  procedure  is  rather  ad  hoc  compared  to  the  wave 
propagation  approach  taken  here,  but  presumably  has  the  advantage  of  being  much  cheaper.  A  major 
shortcoming  of  the  present  method  is  its  expense.  If  tangential  wave  propagation  is  used  in  two  dimensions, 
then  at  each  cell  interface  m  +  1  Riemann  problems  must  be  solved  and  waves  propagated.  Even  using 
Roe’s  approximate  Riemann  solver  (which  is  easy  to  vectorize  over  multiple  Riemann  problems,  particularly  on 
a  Cartesian  grid),  this  may  be  prohibative  for  a  practical  algorithm.  On  the  other  hand,  it  does  not  seem  that  so 
much  information  should  be  necessary,  and  efforts  are  currently  under  way  to  find  simplifications  of  the 
algorithm  by  using  a  smaller  number  of  approximate  waves  that  carry  the  essential  information.  One  possibility 
is  to  use  a  different  kind  of  wave  decomposition  in  place  of  the  Riemann  solutions.  Some  of  the  ideas  of 
Roe[24]  or  Deconinck,  Hirsch  and  Peuteman[8]  might  be  useful  in  this  regard. 

Another  practical  possibility  is  to  use  die  algorithm  developed  here  only  at  irregular  boundaries  of  a 
Cartesian  grid,  coupling  it  to  a  more  efficient  standard  method  in  the  regular  interior  cells. 

On  the  other  hand,  in  some  applications  the  approach  discussed  here  may  have  advantages  over  standard 
finite  volume  or  dimensionally  split  methods  even  on  regular  grids.  The  fact  that  waves  are  propagated  oblique 
to  the  grid  (when  tangential  propagation  is  included)  may  be  advantageous  in  modeling  complex  two 
dimensional  phenomena.  Tests  have  shown  that  including  tangential  propagation,  even  when  not  required  for 
stability,  increases  the  accuracy  with  which  waves  oblique  to  the  grid  can  be  tracked. 

In  some  special  applications  it  could  also  prove  useful  to  use  Courant  numbers  larger  than  one  even  on  a 
regular  grid.  For  example,  in  problems  with  widely  varying  wave  speeds  where  the  phenomena  of  primary 
interest  occurs  on  the  slow  scale,  the  fast  waves  may  carry  little  information  bu*  v-i!!  limit  the  time  step  with 
traditional  methods.  Using  the  wave  propagation  approach  allows  one  to  choose  a  time  step  appropriate  for  the 
slow  scale  phenomena.  This  could  lead  to  a  simultaneous  increase  in  effeiciency  and  reduction  in  the  numerical 
dissipation  of  slow  waves. 
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These  possibilities  are  currently  under  study  and  will  be  reported  on  elsewhere. 

Acknowledgment  It  is  a  pleasure  to  thank  Marsha  Berger  for  numerous  valuable  conversations. 

References. 

[1]  M.  Berger,  thesis,  Stanford  Computer  Science  DepL  Report  No.  STAN-CS-82-924,  1982  (unpublished). 

[2]  M.  Berger,  ICASE  Report  No.  84-43  (1984),  to  appear  in  SINUM. 

[3]  J.P.  Boris  and  DI..  Book,  J.  Comput.  Phys.  11(1973),  38. 

[4]  Y.  Brenier,  SIAMJ.  Numer.  Anal.  21(1984).  1013-1037. 

[5]  I.L.  Chem  and  P.  Colella,  Lawrence  Livermore  National  Laboratory  Report  No.  UCRL-97200  (1987). 

[6]  D.K.  Clarke,  M.D.  Salas,  and  H.A.  Hassan,  AlAA  Journal  24(1986),  353-358. 

[7]  S.F.  Davis.  NASA  Contractor  Report  No.  172179,  ICASE  Report  No.  83-37  (1983). 

[8]  H.  Dcconinck,  C.  Hirsch,  and  J.  Peuteman,  in  10th  Int.  Conf.  on  Num.  Meth.  in  Fluid  Dyn.,  Beijing, 
(1986). 

[9]  H.M.  Glaz,  P.  Colella.  LI.  Glass,  and  R.L.  Deschambault  Proc.  R.  Soc.  bond.  A  398(1985),  117-140. 

[10]  J.B.  Goodman  and  RJ.  LeVeque,  NASA  Contractor  Report  No.  177994,  ICASE  Report  No.  85-44  (1985), 
to  appear  in  SLAM  J.  Numer.  AnaL 

[1 1]  A.  Harten,  J.  Comput.  Phys.  49(1983),  357-393. 

[121  A.  Harten  and  J.M.  Hyman,  J.  Comput.  Phys.  (1983) . 

[13]  P.D.  Lax,  Hyperbolic  systems  of  conservation  laws  and  the  mathematical  theory  of  shock  waves,  SIAM 
Regional  Conference  Series  in  Applied  Mathematics  #11,  (1973). 

[14]  RJ.  LeVeque,  Comm.  Pure  Appl.  Math.  37(1984),  463-477. 

[15]  R  J.  LeVeque,  SIAM  J.  Numer.  Anal.,  22(1985),  1051-1073. 

[16]  RJ.  LeVeque,  in  Proc.  7th  Int.  Corf.  Comp.  Methods  in  Appl.  Sci.  and  Eng.,  Versailles,  France,  (1985). 

[17]  RJ.  LeVeque,  SIAM  J.  Numer.  Anal,  (to  appear). 

[18]  RJ.  LeVeque  and  J.B.  Goodman,  Lectures  m  App/.  Aforh.  22(1985),  51-62. 

[19]  S.  OdJcr  and  S.  Chakravarthy.  SIAM  J.  Numer.  Anal.  21(1984),  955-984. 


-19- 

[20]  P.L.  Roc,  in  Proc.  Seventh  Inti.  Conf.  on  Hum.  Meth.  in  Fluid  Dynamics,  edited  by  W.C.  Reynolds  and 
R.W.  MacCormack,  Lecture  Notes  in  Physics  No.  141,  (Springer,  1981). 

[21]  P.L.  Roe,  /.  Comput.  Phys.  43(1981).  357-372. 

[22]  P.L.  Roc,  in  Numerical  Methods  for  fluid  dynamics,  edited  by  K.W.  Morton  and  M.J.  Baines,  (Academic 
Press,  1982). 

[23]  P.L.  Roe,  Lectures  in  Appl.  Math.  22(1985),  163-193. 

[24]  P.L.  Roe.  NASA  Contractor  Report  No.  172574,  ICASE  Report  No.  85-18  (1985). 

[25]  P.E.  Rubbert,  J£.  Bussolctti,  et.al.,  in  Computational  Mechanics  -  Advances  and  Trends,  AMD  vol.  75, 
edited  by  A.K.  Noor.  49-84.  (1986). 

[26]  P.K.  Swcby,  SIAM  J.  Numer.  Anal.  21(1984),  995-1011. 

[27]  B.  van  Lea,  J.  Comput.  Phys.  32(1979),  101-136. 

[28]  B.  Wedan  and  J.C.  South.  AIAA  Paper  83-1889-CP,  6th  AIAA  Computational  Fluid  Dynamics  Conf., 
(1983). 

[29]  P.  Woodward  and  P.  Colella,  J.  Comput.  Phys.  54(1984),  1 15-173. 

[30]  S.T.  Zalcsak,  J.  Comp.  Phys.  31(1979),  335. 


-22- 


1.10 


1.05 


1.00 


0.95 


0.90 


Figure  5.2.  Density  in  a  smooth  solution  computed  on  a  uniform  grid  (solid 
line)  and  on  nonuniform  grid  (squares)  with  Courant  number  2500.  The  initial  con¬ 
ditions  p(x,0)  are  also  shown. 


Figure  6.1.  Two  grid  cells  Ci  and  C,  from  an  unstructured  gri±  Wave  propaga¬ 
tion  from  the  cell  interface  is  determined  by  solving  a  one  dimensional  Riemann 
problem  in  the  ^-direction,  giving  rise  to  waves  propagating  with  velocities  sj,  sz, 
and  S3. 


figure  6.2.  Propagation  of  a  wave  over  a  general  unstructured  grid.  Cell 
values  are  updated  in  every  cell  the  wave  intersects.  The  contribution  to  each  cell 
is  proportional  to  the  ratio  of  the  shaded  area  within  that  cell  to  the  total  cell  area. 


Figiire  6.3.  Wave  propagation  for  the  scalar  problem  (6.12)  when  tangential 
wave  propagation  is  included.  Shown  for  a  time  step  k  =  h/2. 


Figure  6.4.  The  wave  .R**  is  split  tangentially  into  two  waves  and  Rp<®> 
with  tangential  speeds  Sjf*^  and 


boundary 


Figure  0.1.  Cartesian  mesh  cut  by  a  flat  boundary.  Reflecting  the  interior  cells 
through  this  boundary  gives  the  fictitious  exterior  cells.  Each  wave  propagating 
across  the  boundary  from  an  interior  cell  interface  is  matched  by  a  wave  propagat¬ 
ing  into  the  region  from  a  fictitious  interface.  This  can  also  be  viewed  as  a 
reflection  of  the  interior  wave  at  the  boundary. 


Flcure  9.2  Density  in  the  shock  tube  along  the  lower  wall.  Values  m 
dary  cells  are  plotted  versus  arclength.  The  solid  line  is  the  true  solution, 
order  method,  (b)  With  high  resolution  correcUon  waves  mcluded. 
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